G4 CUT&Tag analysis mESC (WT, FANCJ KO, DHX36 KO, DKO)

Simon Elsässer, Karolinska Institutet (2023)

DESeq2 analysis of peak universe (joint peak set called in any condition, must be called in all three replicates, MACS3)

bw_dir <- "/Volumes/DATA/DATA/Puck/bigwig/"
chromhmm <- "../genome/ESC_10_segments.mm39.bed"

library("wigglescout")
library("ggpubr")
library("ggplot2")
library("DESeq2")
library("dplyr")
library("ggrastr")

clean <- function (fn) {
  fn <- gsub(pattern = ".+/", "", x = fn)
  fn <- gsub(pattern = ".mm9.+", "", x = fn)
  fn <- gsub(pattern = ".mm39.+", "", x = fn)
  fn <- gsub(pattern = "_S.+", "", x = fn)
  fn <- gsub(pattern = ".scaled.bw", "", x = fn)
  fn <- gsub(pattern = ".unscaled.bw", "", x = fn)
  fn <- gsub(pattern = "_batch2", "", x = fn)
  fn <- gsub(pattern = "-", " ", x = fn)
  fn <- gsub(pattern = "_", " ", x = fn)
  fn <- gsub(pattern = " HA ", " ", x = fn)
  fn <- gsub(pattern = "D1D6", "FANCJ-/-", x = fn)
  fn <- gsub(pattern = "P2D2", "DHX36-/-", x = fn)
  fn <- gsub(pattern = "P3D4", "FANCJ-/-DHX36-/-", x = fn)
  return(fn)
}

BWs <- paste0(bw_dir,list.files(bw_dir,pattern="G4_.+R.\\.bw"))

mypal <-c("cornflowerblue","orange","red2","#505050")
mypal3 <-c("cornflowerblue","cornflowerblue","cornflowerblue","orange","orange","orange","red2","red2","red2","505050","505050","505050")

bw_granges_diff_analysis <- function(granges_c1,
                                     granges_c2,
                                     label_c1,
                                     label_c2,
                                     estimate_size_factors = FALSE,
                                     as_granges = FALSE) {

  # Bind first, get numbers after
  names_values <- NULL
  fields <- names(mcols(granges_c1))

  if ("name" %in% fields) {
    names_values <- mcols(granges_c1)[["name"]]
    granges_c1 <- granges_c1[, fields[fields != "name"]]
  }

  fields <- names(mcols(granges_c2))
  if ("name" %in% fields) {
    granges_c2 <- granges_c2[, fields[fields != "name"]]
  }

  cts_df <- cbind(data.frame(granges_c1), mcols(granges_c2))

  if (! is.null(names_values)) {
    rownames(cts_df) <- names_values
  }

  # Needs to drop non-complete cases and match rows
  complete <- complete.cases(cts_df)
  cts_df <- cts_df[complete, ]

  values_df <- cts_df[, 6:ncol(cts_df)] %>% dplyr::select(where(is.numeric))
  cts <- get_nreads_columns(values_df)

  condition_labels <- c(rep(label_c1, length(mcols(granges_c1))),
                        rep(label_c2, length(mcols(granges_c2))))


  coldata <- data.frame(colnames(cts), condition = as.factor(condition_labels))

  dds <- DESeq2::DESeqDataSetFromMatrix(countData = cts,
                                colData = coldata,
                                design = ~ condition,
                                rowRanges = granges_c1[complete, ])


  if (estimate_size_factors == TRUE) {
    dds <- DESeq2::estimateSizeFactors(dds)
  }
  else {
    # Since files are scaled, we do not want to estimate size factors
    sizeFactors(dds) <- c(rep(1, ncol(cts)))
  }

  dds <- DESeq2::estimateDispersions(dds)
  dds <- DESeq2::nbinomWaldTest(dds)

  if (as_granges) {
    result <- DESeq2::results(dds, format = "GRanges",alpha = 0.01)
    if (!is.null(names_values)) {
      result$name <- names_values[complete]
    }

  }
  else {
    # result <- results(dds, format="DataFrame")
    result <- dds
  }

  result
}

get_nreads_columns <- function(df, length_factor = 100) {
  # Convert mean coverages to round integer read numbers
  cts <- as.matrix(df)
  cts <- as.matrix(cts[complete.cases(cts),])
  cts <- round(cts*length_factor)
  cts
}

peak_universe <- "../peaks/G4_combined_min3rep.bed"

BWs.WT <- BWs[grep("WT",BWs)]
BWs.FANCJ <- BWs[grep("D1D6",BWs)]
BWs.DHX36 <- BWs[grep("P2D2",BWs)]
BWs.DKO <- BWs[grep("P3D4",BWs)]
  
# Calculate here some loci or bins
cov.WT <- bw_loci(BWs.WT, loci = peak_universe)
cov.DHX36 <- bw_loci(BWs.DHX36, loci = peak_universe)
cov.FANCJ <- bw_loci(BWs.FANCJ, loci = peak_universe)
cov.DKO <- bw_loci(BWs.DKO, loci = peak_universe)

cov.WT$name <- paste0("peak_",1:length(cov.WT))
cov.DHX36$name <- paste0("peak_",1:length(cov.DHX36))
cov.FANCJ$name <- paste0("peak_",1:length(cov.FANCJ))
cov.DKO$name <- paste0("peak_",1:length(cov.DKO))

diff_DHX36 <- bw_granges_diff_analysis(cov.WT, cov.DHX36,
                                     "WT", "DHX36KO")

converting counts to integer mode
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates

diff_FANCJ <- bw_granges_diff_analysis(cov.WT, cov.FANCJ,
                                     "WT", "FANCJ")

converting counts to integer mode
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates

diff_DKO <- bw_granges_diff_analysis(cov.WT, cov.DKO,
                                     "WT", "DKO")

converting counts to integer mode
gene-wise dispersion estimates
mean-dispersion relationship
-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.
final dispersion estimates

# This takes care of low conts and things like this, but you can also use
# diff_results as is for the things below
lfc_DHX36 <- DESeq2::lfcShrink(diff_DHX36, coef = "condition_WT_vs_DHX36KO", type="apeglm")

using 'apeglm' for LFC shrinkage. If used in published research, please cite:
    Zhu, A., Ibrahim, J.G., Love, M.I. (2018) Heavy-tailed prior distributions for
    sequence count data: removing the noise and preserving large differences.
    Bioinformatics. https://doi.org/10.1093/bioinformatics/bty895

lfc_FANCJ <- DESeq2::lfcShrink(diff_FANCJ, coef = "condition_WT_vs_FANCJ", type="apeglm")

using 'apeglm' for LFC shrinkage. If used in published research, please cite:
    Zhu, A., Ibrahim, J.G., Love, M.I. (2018) Heavy-tailed prior distributions for
    sequence count data: removing the noise and preserving large differences.
    Bioinformatics. https://doi.org/10.1093/bioinformatics/bty895

lfc_DKO <- DESeq2::lfcShrink(diff_DKO, coef = "condition_WT_vs_DKO", type="apeglm")

using 'apeglm' for LFC shrinkage. If used in published research, please cite:
    Zhu, A., Ibrahim, J.G., Love, M.I. (2018) Heavy-tailed prior distributions for
    sequence count data: removing the noise and preserving large differences.
    Bioinformatics. https://doi.org/10.1093/bioinformatics/bty895

data_DHX36 <- plotMA(lfc_DHX36, returnData = T)
data_DHX36$lfc <- -data_DHX36$lfc
data_DHX36$mean <- log10(data_DHX36$mean)

data_FANCJ <- plotMA(lfc_FANCJ, returnData = T)
data_FANCJ$lfc <- -data_FANCJ$lfc
data_FANCJ$mean <- log10(data_FANCJ$mean)

data_DKO <- plotMA(lfc_DKO, returnData = T)
data_DKO$lfc <- -data_DKO$lfc
data_DKO$mean <- log10(data_DKO$mean)

ggscatter(data_DHX36,x ="mean",y="lfc",color="isDE",size = 0.8, alpha=0.5, palette = c("gray",mypal[1])) + geom_hline(yintercept = 0, linetype="dashed", size=0.1) + coord_cartesian(ylim=c(-5,5))

plot_MA_DHX36 <- rasterize(last_plot(), layers='Point', dpi=600) 

ggscatter(data_FANCJ,x ="mean",y="lfc",color="isDE",size = 0.8, alpha=0.5, palette = c("gray",mypal[2])) + geom_hline(yintercept = 0, linetype="dashed", size=0.1) + coord_cartesian(ylim=c(-5,5))

plot_MA_FANCJ <- rasterize(last_plot(), layers='Point', dpi=600) 

ggscatter(data_DKO,x ="mean",y="lfc",color="isDE",size = 0.8, alpha=0.5, palette = c("gray",mypal[3])) + geom_hline(yintercept = 0, linetype="dashed", size=0.1) + coord_cartesian(ylim=c(-5,5))

plot_MA_DKO <- rasterize(last_plot(), layers='Point', dpi=600) 

data_DKO$lfc_DHX36 <- data_DHX36$lfc
data_DKO$lfc_FANCJ <- data_FANCJ$lfc
ggscatter(data_DKO,x ="lfc_DHX36",y="lfc",color="isDE",size = 0.8, alpha=0.5, palette = c("gray",mypal[3]))+ geom_hline(yintercept = 0, linetype="dashed", size=0.1) + geom_vline(xintercept = 0, linetype="dashed", size=0.1) + coord_cartesian(ylim=c(-5,5), xlim=c(-5,5))

plot_LFC_DKO_vs_DHX36 <- rasterize(last_plot(), layers='Point', dpi=600)

ggscatter(data_DKO,x ="lfc_FANCJ",y="lfc",color="isDE",size = 0.8, alpha=0.5, palette = c("gray",mypal[3]))+ geom_hline(yintercept = 0, linetype="dashed", size=0.1) + geom_vline(xintercept = 0, linetype="dashed", size=0.1) + coord_cartesian(ylim=c(-5,5), xlim=c(-5,5))

plot_LFC_DKO_vs_FANCJ <- rasterize(last_plot(), layers='Point', dpi=600)

data_DHX36$lfc_DKO <- data_DKO$lfc
ggscatter(data_DHX36,x ="lfc",y="lfc_DKO",color="isDE",size = 0.8, alpha=0.5, palette = c("gray",mypal[1]))+ geom_hline(yintercept = 0, linetype="dashed", size=0.1) + geom_vline(xintercept = 0, linetype="dashed", size=0.1) + coord_cartesian(ylim=c(-5,5), xlim=c(-5,5))

bed <- as.data.frame(cov.DHX36)[,1:3]

results_table <- data.frame(chr=bed$seqnames,start=bed$start,end=bed$end,name=cov.WT$name, baseMean=lfc_DKO$baseMean,DHX36=data_DHX36$isDE,FANCJ=data_FANCJ$isDE,DKO=data_DKO$isDE,DHX36lfc=data_DHX36$lfc,FANCJlfc=data_FANCJ$lfc,DKOlfc=data_DKO$lfc,dir=".")

results_table$nonsig <- !(results_table$DHX36 | results_table$FANCJ | results_table$DKO)

results_table$sig <- factor("lowlfc", levels=c("nonsig","lowlfc","FANCJ","DHX36","DKO"))
results_table$sig[(results_table$nonsig)] <- "nonsig"
results_table$sig[(results_table$DKO & results_table$DKOlfc>1)] <- "DKO"
results_table$sig[(results_table$DHX36 & results_table$DHX36lfc>1)] <- "DHX36"
results_table$sig[(results_table$FANCJ & results_table$FANCJlfc>1)] <- "FANCJ"

write.table(results_table,"G4_CnT_combined_peaks_DESeq_results.txt", row.names = F, col.names = T,quote = F, sep = "\t")

results_table$sig <- factor(levels=c("nonsig","lowlfc","FANCJ","DHX36","DKO"))

Error in `$<-.data.frame`(`*tmp*`, sig, value = integer(0)) : 
  replacement has 0 rows, data has 71445

library(eulerr)
v <- list(DHX36=DHX36_bed$name,FANCJ=FANCJ_bed$name,DKO=DKO_bed$name)
plot_Venn_all <- plot(euler(v),quantities=T, border="black")
plot_Venn_all

library(eulerr)
DKOneg <- results_table[(results_table$DKO==TRUE & results_table$DKOlfc< -1),c(4)]
v <- list(all=cov.WT$name,DKOneg=DKOneg,DKO=DKO_bed$name)
plot_Venn_DKO <- plot(euler(v),quantities=T)
plot_Venn_DKO

dir.create("./plots",showWarnings = F)
ggsave("plots/peaks_DESeq2_MA_DHX36.pdf",plot_MA_DHX36,width = 4, height= 4)
ggsave("plots/peaks_DESeq2_MA_FANCJ.pdf",plot_MA_FANCJ,width = 4, height= 4)
ggsave("plots/peaks_DESeq2_MA_DKO.pdf",plot_MA_DKO,width = 4, height= 4)
ggsave("plots/peaks_DESeq2_LFC_DHX36.pdf",plot_LFC_DKO_vs_DHX36,width = 4, height= 4)
ggsave("plots/peaks_DESeq2_LFC_FANCJ.pdf",plot_LFC_DKO_vs_FANCJ,width = 4, height= 4)
ggsave("plots/peaks_DESeq2_Venn_KOs.pdf",plot_Venn_all,width = 4, height= 4)
ggsave("plots/peaks_DESeq2_Venn_DKO_vs_all.pdf",plot_Venn_DKO,width = 4, height= 4)

library(cowplot)
dir.create("./panels",showWarnings = F)

library("cowplot")
p <- ggdraw() +
  draw_plot(plot_MA_DHX36, x = 0, y = 0.5, width = .33, height = .5) +
  draw_plot(plot_MA_FANCJ, x = .33, y = 0.5, width = .33, height = .5) +
  draw_plot(plot_MA_DKO, x = 0.66, y = 0.5, width = .33, height = 0.5) +
  draw_plot(plot_LFC_DKO_vs_DHX36, x = 0, y = 0, width = 0.33, height = 0.5) +
  draw_plot(plot_LFC_DKO_vs_FANCJ, x = 0.33, y = 0, width = 0.33, height = 0.5) +
  draw_plot(plot_Venn_all, x = 0.66, y = 0, width = 0.33, height = 0.5) +
  draw_plot_label(label = c("a", "b", "c","d","e","f"), size = 15,
                  x = c(0, 0.33, 0.66, 0, 0.33, 0.66), y = c(1, 1, 1, 0.5, 0.5, 0.5))
p

ggsave("panels/peaks_DESeq2.pdf",p,width=10, height=6)

library(reshape2)
mdf <- melt(data.frame(name=rownames(cov),cov[,6:21]))

Using name, DHX36_sig, FANCJ_sig, DKO_sig, non_sig as id variables

mdf <- mdf[mdf$value<500,]
mdf$cond <- "WT"
mdf$cond[grep("DHX",mdf$variable)] <- "DHX36-/-"
mdf$cond[grep("FAN",mdf$variable)] <- "FANCJ-/-"
mdf$cond[grep("DKO",mdf$variable)] <- "DHX36-/-FANCJ-/-"
plot_viol_rep_all_peaks <- ggviolin(mdf, x="variable",y="value",fill="cond",palette = mypal, add="mean_sd") +
  coord_cartesian(ylim=c(0,30))
plot_viol_rep_all_peaks

plot_viol_rep_DHX36_peaks <- ggviolin(mdf[mdf$DHX36_sig,], x="variable",y="value",fill="cond",palette = mypal, add="mean_sd") +
  coord_cartesian(ylim=c(0,30))
plot_viol_rep_DHX36_peaks

ggviolin(mdf[mdf$FANCJ_sig,], x="variable",y="value",fill="cond",palette = mypal, add="mean_sd") +
  coord_cartesian(ylim=c(0,30))

plot_viol_rep_DKO_peaks <- ggviolin(mdf[mdf$DKO_sig,], x="variable",y="value",fill="cond",palette = mypal, add="mean_sd") +
  coord_cartesian(ylim=c(0,30))
plot_viol_rep_DKO_peaks

plot_viol_rep_nonsig_peaks <- ggviolin(mdf[mdf$non_sig,], x="variable",y="value",fill="cond",palette = mypal, add="mean_sd") +
  coord_cartesian(ylim=c(0,30))
plot_viol_rep_nonsig_peaks

ggsave("plots/peaks_DESeq2_viol_rep_all.pdf",plot_viol_rep_all_peaks,width = 2, height= 3)
ggsave("plots/peaks_DESeq2_viol_rep_DKO.pdf",plot_viol_rep_DKO_peaks,width = 2, height= 3)
ggsave("plots/peaks_DESeq2_viol_rep_DHX36.pdf",plot_viol_rep_DHX36_peaks,width = 2, height= 3)
ggsave("plots/peaks_DESeq2_viol_rep_nonsig.pdf",plot_viol_rep_nonsig_peaks,width = 2, height= 3)

sdf <- aggregate(value ~ name + cond, data=mdf, FUN="mean")
sdf$cond <- factor(sdf$cond,levels=c("WT","DHX36-/-","FANCJ-/-","DHX36-/-FANCJ-/-"))
ggviolin(sdf, x="cond",y="value",fill="cond",palette = mypal[c(4,1,2,3)], add="mean_sd") +
  coord_cartesian(ylim=c(0,50))

sdf$DHX36_sig <- sdf$name %in% DHX36_bed$name
sdf$FANCJ_sig <- sdf$name %in% FANCJ_bed$name
sdf$DKO_sig <- sdf$name %in% DKO_bed$name
sdf$DKO_top <- sdf$name %in% DKO_bed_top$name
ggviolin(sdf[sdf$DHX36,], x="cond",y="value",fill="cond",palette = mypal[c(4,1,2,3)], add="mean_sd") +
  coord_cartesian(ylim=c(0,75))

ggviolin(sdf[sdf$FANCJ_sig,], x="cond",y="value",fill="cond",palette = mypal[c(4,1,2,3)], add="mean_sd") +
  coord_cartesian(ylim=c(0,75))

ggviolin(sdf[sdf$DKO_sig,], x="cond",y="value",fill="cond",palette = mypal[c(4,1,2,3)], add="mean_sd") +
  coord_cartesian(ylim=c(0,75))

ggviolin(sdf[sdf$DKO_top,], x="cond",y="value",fill="cond",palette = mypal[c(4,1,2,3)], add="mean_sd") +
  coord_cartesian(ylim=c(0,75))

annotate peak table

write.table(results_table,"G4_CnT_combined_peaks_DESeq_annotated.txt", row.names = F, col.names = T,quote = F, sep = "\t")

Error: unexpected symbol in:
"
write.table"

#read in - if you want to skip de novo generation
bed_peaks_table <- read.table("G4_CnT_combined_peaks_DESeq_annotated.txt", header = T, sep = "\t")

vpal=colorRampPalette(c("cornflowerblue","orange","red2"))
stats <- as.data.frame(table(bed_peaks_table$V4))
ggdonutchart(stats,x = "Freq",label="Var1",fill=vpal(10))

p <- ggdraw() +
  draw_plot(plot_bar_chromhmm_anno, x = 0, y = 0, width = .20, height = 1) +
  draw_plot(plot_viol_rep_DKO_peaks, x = .2, y = 0, width = .40, height = 0.7) +
  draw_plot(plot_viol_rep_nonsig_peaks, x = 0.6, y = 0, width = .40, height = 0.7) +
  draw_plot_label(label = c("a", "b", "c"), size = 15,
                  x = c(0, 0.2, 0.6), y = c(1, 0.8, 0.8))
p

ggsave("panels/peak_anno_violin_v1.pdf",p,width=8, height=4)

p <- ggdraw() +
  draw_plot(plot_bar_chromhmm_anno, x = 0, y = 0, width = .4, height = 1) +
  draw_plot(plot_viol_rep_DKO_peaks, x = .4, y = 0.5, width = .60, height = 0.5) +
  draw_plot(plot_viol_rep_nonsig_peaks, x = 0.4, y = 0, width = .60, height = 0.5) +
  draw_plot_label(label = c("a", "b", "c"), size = 15,
                  x = c(0, 0.4, 0.4), y = c(1, 1, 0.5))
p

ggsave("panels/peak_anno_violin_v2.pdf",p,width=8, height=3)